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StCTION  1 
INTRODUCTION 


Waves  are  influenced  by  the  constitutive  properties  of  the  material 
in  which  they  propagate.  The  effects  of  constitutive  properties  on  spherical 
wave  propagation  is  of  special  interest  in  calculating  the  effects  of  under¬ 
ground  nuclear  explosions.  Recent  efforts  using  numerical  methods  to  predict 
the  effects  of  the  PILEDRIVER  event  have  focussed  attention  on  the  role  of 
mathematical  models  of  constitutive  properties.  The  purpose  of  the  present 
paper  is  to  illustrate  how  various  assumptions  about  constitutive  properties 
affect  the  calculated  wave  propagation. 

In  the  present  analysis,  the  finite  element  method  is  adapted  to 
spherical  geometry.  A  mathematical  model  is  derived  for  Cedar  City  Tonal ite 
based  on  laboratory  measurements.  The  model  contains  a  bulk  modulus  which 
depends  on  the  current  density,  a  shear  modulus,  a  yield  criterion  and  a  rule 
of  plastic  flow.  The  Tonalite  is  considered  to  be  inf'nite  medium  surround¬ 
ing  a  cavity  which  contains  a  sphere  of  chemical  explosive.  Detonation  of  the 
explosive  is  represented  by  applying  a  pressure  to  the  surface  of  the  cavity 
which  varies  with  time  in  a  manner  similar  to  that  measured  by  Physics  Inter¬ 
national  Company  (Reference  1).  The  peak  pressure  used  in  the  present  calcu¬ 
lations  is  31 .5  kilobars. 

The  results  of  the  calculations  are  presented  as  stress/time  histories 
and  stress/strain  relations  at  various  ranges  and  as  rates  of  attenuation  of 
the  peak  radial  stress.  The  effects  of  varying  these  properties  are  studied 
by  comparing  stress/time  histories  and  stress/strain  relations  at  various 
ranges  and  as  rates  of  attenuation  of  the  peak  radial  stress  The  primary 
comparisons  are  made  among  cases  where  the  radial  stress  applied  to  the  cavity 
surface  and  hence  the  impulse  is  invariant. 
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FINITE  ELEMENT  METHOD 


The  application  of  the  finite  element  method  to  problems  in 
continuum  mechanics  has  been  thoroughly  discussed  by  previous  writers 
(References  2-4).  Hence,  only  a  few  subjects  pertaining  specifically  to  the 
present  adaptation  of  the  finite  element  method  are  discussed  below. 

The  incremental  equation  of  dynamic  equilibrium  is 


[m] { 6u }  +  [c] { fiu)  +  [K]{6u}  -  {6P}  0) 

*  Global  mass  matrix 

-  Global  damping  matrix 
■  Global  (tangent)  stiffness  matrix 

-  Increments  of  acceleration,  velocity,  and  displacement 
■=  Increment  of  load 

The  global  stiffness  matrix  is  assembled  from  element  stiffness 
matrices  [k] .  For  a  typical  element  (see  Figure  2-1) 


in  which  [m] 
[c] 
[K] 

{6u}.  { 6u} ,  {6u} 
{6P} 


[k] 


[B]T[D] [B] dV [A] ” 1 


(2) 


% 


where  [A]  '  =  A  transformation  matrix  relating  the  generalized  displacements 

of  the  finite  element  theory  to  the  physical  radial  displace¬ 
ment 


[B]  =  Strain/displacement  transformation  matrix 

[D]  =  Matrix  of  tangent  stress/strain  moduli  relating  incremental 

stress  to  incremental  strain 


PRECEDE  M6E  BLANK 
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FIGURE  2-1.  SPHERICAL  ELEMENT  DEFINED  BY 
TH  TH 

I  AND  J  NODAL  SURFACES 


In  the  present  case 


•  « 

_L  _L 

2  '2 

,  ri  r' 

1l.1l 

?.  2 

r.  r.  -r,  r. 

J  1  L  j  i  J 


where  rj*  rj  “  Radii  of  the  i  ^  and  nodal  surfaces 

[A]  1  expresses  the  displacement  u  between  the  Ith  and  jth  nodal 
surfaces  by  the  following  equation: 


u  ■  <!r  r"2>  [A]’1  j  1 

(  j 


The  displacement  function.  Equation  2-k,  satisfies  the  elasto-stat ic 
solution  of  a  hollow  sphere  subjected  to  internal  and/or  external  pressure. 
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FIGURE  2-2.  LUMPED  MASS 
APPROXIMATION 
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steps  as  illustrated  in  Figure  2-3.  In  the  first  step,  the  solution  at 

t  +  6t/2  is  estimated  using  the  known  tangent  stiffness  at  t  .  This  is  done 
o  o 

by  obtaining  a  temporary  solution  at  t  +  6t  (1  in  the  Figure)  and  then  inter¬ 
polating  to  tQ  +  6t/2  (2  in  the  Figure).  The  solution  at  the  end  of  the  first 

pass  is  used  to  evaluate  the  stiffness  matrix  [K]  at  t  +  6t/2.  The  first 

o 

pass  solution  is  then  discarded.  The  second  pass  is  performed  using  [K]  at 
tQ  +  6t/2  to  obtain  a  temporary  solution  at  tQ  +  26t  (3  in  the  Figure).  The 

final  result  is  obtained  by  interpolating  between  t  and  t.Q  +  26t  to  find  the 
incremental  displacement  at  t  +  fit  ( k  in  the  Figure). 


aaritss 


naroRMtv 

sainiw - 
M  PASS 


FIGURE  2-3.  INTEGRATION  TECHNIQUE 
FOR  PRESENT  FINITE  ELEMENT  METHOD 


This  integration  procedure  requires  the  equations  of  motion  to  be 
solved  twice  during  every  time  step  and  therefore  appears  to  be  more  time- 
consuming  than  necessary.  However,  the  materials  of  interest  exhibit  a  large 
difference  in  stiffness  between  loading  and  unloading,  and  it  was  feared  that 
small  numerical  errors  might  lead  to  large  changes  in  stiffness  and  to  insta¬ 
bility.  Hence,  the  conservative  procedure  described  above  was  adopted.  The 
present  calculations  show  that  this  procedure  is  stable.  However,  some  recent 
theoretical  work  at  AJA  seems  to  indicate  that  a  one  pass  method  would  be 
just  as  effective  as  the  present  two  pass  method.  Further  work  is  required 
to  conf i rm  this . 
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MATHEMATICAL  MODEL  OF  MATERIAL  PROPERTIES 


The  mathematical  model  which  is  used  in  the  present  work  has  the 
following  basic  features: 

a.  Bulk  modulus  (B)  may  be  a  function  of  the  excess  compression 

p 

V  ■  p  -  1,  where  p  ■  current  density  and  PQ  •  initial 

o 

dens i ty . 

b.  Shear  modulus  (6)  may  be  a  function  of  the  current  state  of 
stress.  In  the  present  examples  it  is  assumed  to  be  constant. 

c.  Yield  criterion  may  be  a  function  of  the  first  stress  invariant 
(J I )  and  of  the  second  invariant  of  stress  deviator  (J£) 
(Reference  6) . 

d.  Work-hardening  or  strain-hardening  rules  prescribe  how  the 
yield  criterion  may  vary  as  a  function  of  plastic  work  or 
plastic  strain.  In  the  present  examples,  hardening  is  assumed 
to  be  zero  and  the  initial  yield  criterion  is  a  permanent 
property  of  the  material  (Reference  6). 

e.  Flow  rules  prescribe  how  changes  in  plastic  strain  are  related 
to  changes  in  stress  when  the  yield  criterion  is  satisfied. 

Two  flow  rules,  the  plastic  potential  (Reference  7)  and  a  vers¬ 
ion  of  the  Prandt 1 -Reuss  rule  (Reference  8),  are  considered. 


This  mathematical  model  must  be  expressed  as  a  matrix  of  coefficients 
| D ]  relating  stress  increments  jdoj  to  total  strain  increments  jdej  . 


The  |D|  matrix  has  two  purposes  in  the  finite  element  scheme  described  above. 
One  is  in  formulating  the  element  stiffness  matrix  [k] ,  Equation  2.  The 
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second  purpose  is  in  determining  the  stress  increments  in  the  element  between 
the  Ith  and  jth  node  from  the  Incremental  nodal  displacements  as  Indicated 
by  Equation  9* 


(do) 


[0]  [B]  [A] 


jdu, 

H 


(9) 


The  starting  point  In  deriving  [0]  Is  Hooke's  law  for  small  strains 
which  relates  ‘he  Incremental  stress  tensor,  do..,  to  the  elastic  component 
of  the  strain  tensor,  de^.. 


do 


ij  ■  id4su  +  2GKj) 


where  X  -  Lame's  parameter,  B  -  *•  G 


(10) 


If  the  state  of  stress  does  not  satisfy  the  yield  criterion,  the  elastic 
strain  increment  and  the  total  strain  Increment  are  equal  and  [D]  Is  based 
on  Equation  10.  However,  If  the  yield  criterion  is  satisfied,  further  com¬ 
putations  are  necessary  to  obtain  [Dj. 


Defining  the  elastic  strain  Increment  as  the  difference  between  the 

p 

plastic  strain  increment,  de{j,  and  the  total  strain  increment,  d e.y 
Equation  10  may  be  rewritten  as  follows: 


do 


Ij 


* (d£kk 


def,  )  6..  +  2G/de.  . 
kk/  ij  \  ij 


(11) 


The  flew  rule 
components  of 


Is  used  to  express  dcjj  in  terms  of  the  yield  criterion  and 
the  stress  or  stress  deviator  tensor 


.  p 
dEij 

A  9f 

(plastic  potential) 

(12a) 

de? , 

•J 

■  *#- 

90  lj 

(present  version  of 
Prandtl-Reuss  flow  rule) 

(12b) 
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For  the  sake  of  simplicity  in  the  following  derivation,  both  Equations  12a,  b 
are  expressed  as 


03) 


In  the  subsequent  applications  of  the  mathematical  model  it  is  made  clear  which 
flow  rule  is  being  used. 


From  the  assumption  of  no  work  or  strain  hardening,  it  follows  that 
there  is  no  change  in  the  yieid  function.  The  mathematical  statement  of  this 
assumption  is 


d (f  (o. j))  *  0 


(14a) 


or 


fi-do.j  =  f..  do..  -  0  (Hb) 

Substituting  Equation  (11)  into  Equation  (14b)  and  making  use  of  Equation  (13) 
leads  to  the  following  equation. 

i(dckk-  Afkk>  VlJ  +  2G<deij  -  Afij>  fij  ■  0  05) 

The  scalar  quantity  A  may  be  found  by  rearranging  Equation  15. 

.  *<d‘kk><fU>  +2G(dc  )(f,) 

Making  use  of  Equation  (16)  in  Equation  (13)  and  substituting  the  result  into 
Equation  (11),  the  stress  increment  is  expressed  in  terms  of  the  total  strain 
increment  and  the  total  stresses.  This  is  the  desired  result. 

The  remaining  task  is  to  factor  out  coefficients  in  the  [ D]  matrix, 
which  is  given  on  the  following  page. 
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where 


Derivatives  of  the  yield  function  f  with  respect 

to  stress  components  (plastic  potential  flow  rule)  or 

stress  deviator  components  (present  version  of 

Prandtl-Reuss  flow  rule),  f  indicates  differentiation 

r 

with  respect  to  radial  stress  or  stress  deviator  component. 
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PROPERTIES  OF  ROCK  USED  IN  FINITE  ELEHENT  CALCULATIONS 

Specific  material  properties  were  derived  from  a  variety  of 
laboratory  tests  on  Cedar  City  Tonal Ite  and  NTS  granite.  Properties  of  both 
rocks  are  Incorporated  into  the  model,  which  is  intended  to  represent  a 
weathered  granite. 

BULK  HODULUS 

The  mathematical  model  of  bulk  modulus  is  based  on  hydrostatic 

compression  tests  In  the  range  0-37  kilobars  (Reference  9,  10).  Release 

adlabat  data  (Reference  11)  were  also  taken  Into  account  when  deciding  that 

hysteresis  should  be  incorporated  Into  the  model.  The  model,  which  is 

intended  for  ur.e  only  when  u  ■  - 1  Is  positive,  is  given  below. 

po 

Loading  (u  >  u  ,  the  previous  maximum  u  In  a  particular  element) 

—  max 


8  ■  "ui,  -  <8ui,  -  V  (^) 


Unloading  or  reloading  (u  >  u  ) 

max 


-  B  +  (B  -  B  )/H-\ 
u  uit  u  I  w2 / 


and  B  ■  the  lesser  of 


I 


B0  +  <Bult  -  V  “f 


(18) 


(19) 


(20) 


HUGONIOT 

HYDROSTAT 


»  REF  14 
•  REF  11 


PRESSURE,  Kb 

I.  LOW  PRESSURE  BUEK  MODULUS  COMPARED  WITH  DATA  IREF  10) 


b  MODEL  HUGONIOT  AND  HYDROSTAT  COMPARED  WITH  DATA 


FIGURE  A-2.  MODEL  BULK  MODULUS,  HYDROSTAT  AND  HUGONIOT  COMPARED  WITH  DATA 
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SHEAR  MODULUS 
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The  model  shear  modulus  Is  based  on  the  slope  of  deviatoric  s cress/ 
strain  curves  measured  during  trlaxla)  compression  tests  on  Cedar  City 
Tonal ite  (Reference  10).  The  data,  which  are  shown  in  Figure  4-3  seem  to 
indicate  that  G  decreases  with  increasing  confining  pressure.  The  trend  is 
too  week  to  Justify  a  model  which  is  more  complicated  then  G  »  a  constant. 
However,  the  G  which  is  chosen 

G  ■  0.69  kb  *  10^  p*i 

gives  more  weight  to  the  data  at  the  high  pressure  end  of  the  range  of 
measurements. 


YIELD  CRITERIA 


EACH  TYPE  OF  SYMBOL  DENOTES  ONE 
TRIAXIAL  COMPRESS  ION  TEST  (REF  10) 


a  a  ?  • 

■  ■  *  1  . 

.♦  ♦■♦A 


PRESENT  MODEL 

_ _ _ _ _ I _ 

0.2  0.4  0.6  0.8  1.0 

PRESSURE,  Kb 

FIGURE  4-3.  SHEAR  MODULUS  VERSUS  PRESSURE 


The  yield  criterion,  which  has  this  form 


■ 


aJ^  -  c  <_  0 


Is  based  on  several  types  of  experimental  data.  The  yield  criterion  at  low 
pressure  is  based  on  static,  triaxlal  compression  experiments  which  were 
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conducted  on  samples  whose  minimum  dimension  was  1  to  2  in.  Data  obtained 
from  presawn  samples  of  Nevada  Test  Site  granite  (Reference  12)  and 
Cedar  City  Tonal ite  (Reference  10)  and  from  intact  samples  of  Westerly 
granite  (Reference  13)  and  Cedar  City  Tonal Ite  (Reference  10)  were  considered. 
The  reason  for  considering  presawn  samples,  which  were  sawn  through  at  angles 
of  45-60  deg  to  the  direction  of  major  principal  stress.  Is  to  try  to  take 
into  account:  preexisting  cracks  ir  the  Physics  International  specimens  and 
the  faults  and  joint  planes  which  exist  in  large  rock  masses  in  situ. 

The  data  on  which  the  yield  criterion  Is  based  are  shown  In 
Figure  4-4. 


FIGURE  4-4.  TRI AXIAL  COMPRESSION  DATA  AND  PRESENT  YIELD  CRITERION 
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A  lower  limit  on  shear  strength  In  the  range  of  the  data  Is  used 
for  the  yield  criterion  in  the  present  calculations.  The  purpose  of  this  is 
to  use  a  physically  plausible  yield  criterion  while  allowing  a  maximum  amount 
of  material  to  undergo  inelastic  deformation.  The  coefficients  defining  the 
present  criterion  are  as  follows. 

aj  ■  -0.167 
Cj  ■  1450  psl 

This  criterion  Is  assumed  to  apply  in  the  following  range  of  mean  stress. 

J, 

-149,350  psl  1  y-  -  p 

The  criterion  implies  that  the  maximum  allowable  deviatoric  stress  Is  zero 
J1 

when  p  -  »  1900  psl.  The  material  has  a  strength  under  uniaxial  tensile 

stress  of  2500  psl,  which  is  too  large  to  be  a  good  respresentatlon  of 
weathered  granite.  A  preliminary  study  indicated  that  the  attenuation  of 
peak  stress  and  general  shape  of  the  compressive  phase  of  the  pulse  varied 
little  as  the  details  of  the  tensile  properties  were  varied.  Hence,  the 
present  Coulomb  yield  criterion  plus  the  restriction  that  the  mean  tensile 
stress  may  never  exceed  1900  psl  were  adopted  for  simplicity. 

The  yield  criterion  which  is  assumed  to  appiv  In  the  region 
Jl  . 

—  ■  p  <  -149,350  psl  is  the  von  Mises  type 


c2  ■  76,000  psl 
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This  criterion  is  based  on  the  hypothesis  that  the  rock  has  an  ultimate  shear 
strength  depending  on  the  strength  of  minerals  withir.  grains  and  is  Independ¬ 
ent  of  friction  between  grain  boundaries.  The  Idea  of  using  a  von  Mines 
yield  criterion  is  further  supported  by  the  data  In  figure  4-2,  whic.i 
Indicate  that  the  hugoniot  is  approximately  a  constant  distance  equai  to 
about  90,000  psl  above  the  hydrostat  for  u  >_0.05.  This  is  consistent  with 
a  value  of  c2  *  76*000  psi. 

The  mean  stress  at  which  the  transition  from  Coulomb  to  von  Mises 
yield  criterion  is  made  is  determined  by  using  the  criterion  giving  the 
minimum  allowable  Vj' 


V-,  ♦  ci 

c2 


One  feature  of  the  model  which  does  not  agree  with  the  experimental 
data  is  the  hugoniot  elastic  limit  (HEl  in  Figure  4-2),  defined  as  the  stress 
level  at  which  a  wave  propagating  in  the  material  has  more  than  one  character¬ 
istic  speed.  Although  this  definition  is  somewhat  inadequate  in  that  it  over¬ 
looks  diffusion,  rate  effects  and  compaction  at  iow  stress  levels,  there  seems 
to  be  a  distinct  stress  level  in  Tonal ite  where  the  wave  spaed  alters  appreciably. 
This  has  been  reported  experimentally  at  about  y  -  0.048  or  Oj  ■  -320,000  psi 
(Reference  14).  In  the  present  model,  it  occurs  at  about  y  ■  0.066  and 
Oj  ■  -435,000  p.ti.  One  way  to  bring  this  aspect  of  the  model  into  better 
agreement  with  the  experimental  data  is  to  increase  the  shear  modulus  to 
G  ■  1.5  x  10'-  psi.  This  would  also  bring  the  shear  modulus  into  better 
agreement  with  low  pressure  data,  and  in  calculations  it  wouid  indicate  more 
inelastic  deformation  on  loading  than  does  the  present  model.  The  overall 
effect  on  the  present  series  of  calculations  of  making  such  a  change  would 
probably  be  small,  however. 
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DENSITY 

It  is  assumed  that  the  initial  density  of  the  material 
•  2.55  gm/cm3  ■  0.000238  It -sec^/in.^.  This  value  is  representative  of  a 
number  of  samples  taken  from  the  Cedar  City  site. 

PREVIOUS  WORK 

The  mathematical  model  of  weathered  granite  used  here  is  adapted 
from  previous  work  reported  in  Reference  15. 
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NUMERICAL  RESULTS 

Calculations  were  performed  In  which  the  cavity  size  and  pressure/ 
time  history  applied  to  the  surface  of  the  cavity  are  comparable  to  those  in 
physical  experiments  conducted  by  Physics  International  Company,  Reference  1. 

The  finite-element  grid  used  in  the  present  calculations  is  shown  in  Figure  5-1a. 
The  pressure/time  history  for  the  present  calculation.  Figure  5“1»  is  adapted 
from  a  radial-stress/time  history  which  was  measured  in  a  spherical  wave 
experiment  on  Cedar  City  Tonal ite  where  the  wave  was  generated  by  detonating 
chemical  explosive. 

The  purpose  of  using  the  fine  mesh  illustrated  in  Figure  5“1a  is  to 

try  to  maintain  the  correct  rise  time  of  the  wave.  The  rise  time,  as  indicated 

by  measurements  in  the  physical  experiment  which  are  reliable  to  the  nearest 

0.1  x  10  6  sec  is  about  0.7  x  10  ^  sec.  As  should  be  expected  in  a  numerical 

solution  method  of  the  type  used  here,  the  rise  time  increases  to  2.0 

*  3.0  x  10  ^  sec  by  the  time  the  wave  has  propagated  a  distance  equal  to 

5  x  r  ,  where  r  is  the  cavity  radius.  A  somewhat  greater  increase  in  rise 

o  o 

time  is  noticed  in  the  physical  experiment,  which  is  probably  due  either  to 
viscous  properties  of  the  material  or  to  dispersion  caused  by  small-scale 
inhomogeneities  in  the  rock  structure. 

The  reason  for  trying  to  represent  the  correct  rise  time  is  that, 

in  the  vicinity  of  the  cavity,  the  rise  time  influences  the  rate  of  attenuation 

of  peak  radial  stress  and  particle  velocity.  This  point  is  discussed  in 

detail  in  Appendix  A.  The  rate  of  attenuation  is  a  parameter  of  practical 

interest  and  is  used  in  the  present  analysis  as  one  index  of  the  effect  of 

changing  a  particular  material  property.  Hence  it  is  important  for  the  rise 

time  in  the  numerical  calculations  to  be  a  physically  meaningful  quantity  and 

for  i t  to  be  the  same  from  one  computation  to  the  next.  Hence,  a  conservative 

-8 

integration  time  step  fit  ■  5  x  10  sec  and  a  fine  finite-element  mesh  were 
used. 

PRECEDING  PAGE  BUNK 
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b.  PRESSURES IME  HISTORY  AT  CAVITY  WALL 


FIGURE  5-1.  MESH  SIZE  AND  LOADING  USED  IN  PRESENT  CALCULATIONS 
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The  main  results  of  the  present  study  were  obtained  from  a  series  of 
seven  c&lculdtions  which  are  described  in  the  foi lowing  table. 


Case 

1 


Description 

Linear  elastic  (no  yield  condition  imposed); 

B  -  Bq  -  1.2  x  106  psl;  G  +  1  x  106  psl 

Variable  modulus  (no  yield  condition  imposed); 
B  defined  by  Equations  18-20;  G  ■  1  x  10^  psi 

Variable  modulus  as  in  Case  2  with  von  Hises 
yield  criterion  (a  *  0,  c  *  76125  psi) 

Variable  modulus  as  in  Case  2  with  Coulomb 
yield  criterion  (o  ■  -0.167,  c  «  1450  psi)  at 
p  >  -l1*. 5  kb  and  von  Mises  criterion  (o  *  0, 
c  ■  76125  psi)  at  p  <  -14.5  kb.  Plastic 
Potential  flow  rule  is  used. 

Same  as  Case  4,  except  that  Prandtl-Reuss 
flow  rule  is  used. 


Same  as  Case  4,  except  that  impulse  is  increased 
by  a  factor  of  2. 

Same  as  Case  2,  except  that  impulse  is  increased 
by  a  factor  of  2. 


STRESS/TIME  HISTORIES 

The  radial  and  circumferential-stress/time  histories  at  various  ranges 
from  the  cavity  wall  are  shown  for  each  case  in  Figures  5”2  to  5"5.  Compressive 
stresses  are  defined  to  be  negative.  The  time  duration  of  each  calculation  is 
35  x  10  ^sec,  which  allows  a  fairly  complete  picture  of  the  stress  pulse  in  the 
range  r/rQ  ■  1-4. 
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FIGURE  5-2.  STRESS/TIME  HISTORIES,  CASES  1  AND  2 
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a.  RADIAL  STRESS 


FIGURE  5-3.  STRESS/TIME  HISTORIES,  CASES  3  AND  A 
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FIGURE  5-4.  STRESS/TIHE  HISTORIES,  CASE  5 
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FIGURE  5-5.  STRESS/TIME  HISTORIES,  CASES  6  AND  7 
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HOOP  TENSILE  STRESS 

The  hoop  stresses  (o0,  o^)  are  tensile  or  compressive  according  to 
the  sense  of  the  elastic  component  of  hoop  strain.  In  spherical  waves,  the 
hoop  strains  (eQ,  e^)  tend  to  become  tensile  due  to  outward  radial  disp'acement 
(u)  according  to  the  following  equation 

e6  “  e*  "  u/r 

Since  in  the  present  calculations  the  rise  time  of  radial  stress  and  particle 
velocity  is  very  short,  u-0  at  the  time  of  peak  The  state  of  strain 

at  this  time  is  approximately  u.iiaxiai  compression  in  the  r-di rection,  and 
hence  compressive  o0,  0^  are  induced  by  Poisson's  ratio  effects  just  as  they 
are  in  plane  wave  propagation.  This  behavior  is  governed  by  the  material 
properties  and  wave  shape  used  in  the  present  calculations.  The  hoop  stress 
would  be  initially  tensile  if  Poisson's  ratio  were  assumed  to  be  zero,  and  the 
amplitude  and  duration  of  the  compressive  hoop  stress  would  be  reduced  if  the 
rise  time  were  lengthened. 

Following  an  initial  compressive  phase,  oQ,  0^  become  tensile  in 
the  linear  elastic  and  variable  modulus  cases  (Cases  1  and  2).  In  Case  3, 
a0,  0^  at  r/r^  -  1.015  are  influenced  by  the  development  of  plastic  deforma¬ 
tion,  during  which  the  effective  Poisson's  ratio  is  0.5.  The  component  of  hoop 
strain  due  to  Poisson's  ratio  effects  is  thus  larger  than  in  Cases  1  and  2,  and 

outward  displacement  must  be  correspondingly  larger  before  a.,  a,  can  become 

n  9 

tensile.  This  accounts  for  the  occurrance  of  tensile  stresses  at  about  t  « 

22  x  10  ^  sec  in  Case  3*  At  ranges  greater  than  about  r/rQ  ■  1.75,  the  yield 

criterion  is  no  longer  satisfied  in  Case  3  and  a,,  0,  become  similar  to  those 

v  9 

in  Cases  1  and  2. 

In  Case  4,  tensile  stresses  do  not  occur  anywhere  within  the  range 
considered.  Part  of  the  reason  for  this  behavior  is  the  Coulomb  yield  condi¬ 
tion  which  allows  plastic  deformations  to  occur  at  relatively  iow  stress  levels 
and  large  ranges.  Also,  dilatency  or  plastic  volume  increase  takes  place  during 
unloading.  In  the  range  of  this  problem,  outward  displacement  is  not  large 
enough  to  absorb  all  this  dilatency,  and  so  hoop  compressive  stresses  are  induced. 
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This  effect  of  di latency  is  further  illustrated  In  Case  5,  where 
the  Coulomb  yield  criterion  is  used  but  plastic  di latency  is  not  allowed. 

The  low  yield  criterion  suppresses  tensile  stresses  at  r/r^  -  1.015,  but 
at  greater  ranges,  tensile  stresses  develop  for  the  same  reason  as  in  Cases  1*-3. 

0/0.  TRAJECTORIES 

r  0 

Plots  of  0_  vs  o.,  o.  further  illustrate  the  influence  of  the 
r  u  9 

yield  criterion  on  the  wave  shape.  Figures  5-6  to  5“8  show  o  /o  ,o 

r  0  ♦ 

trajectories  at  r/rQ  ■  1.015,  2.00  and  3*00.  Superposed  on  each  figure  are 
yield  criteria  used  in  Case  3  (von  Hises)  and  Cases  4,  5  (von  Kises  at 
p  <  I4.5kb;  Coulomb  elsewhere).  In  Cases  1  and  2,  the  stress  path  is  unre¬ 
strained  by  a  yield  criterion.  The  ratio  of  0  /o.  is  greater  in  Case  1  than 

r  d 

in  Case  2  because  Poisson's  ratio  in  Case  1  is  smaller.  In  Case  2,  Poisson's 
ratio  tends  toward  0.45  as  the  bulk  modulus  B  tends  toward  its  maximum  value  of 
7.6  x  106ps i .  This  causes  the  material  in  Case  2  to  appear  more  fluid-iike 
than  in  Case  1.  Hence  tiie  stress  state  in  Case  2  lies  nearer  the  hydrostat. 


FIGURE  5-6.  STRESS  TRAJECTORIES  AT 

r/r  -  1.015 
o 


26 


B5BO 


R-681 3-777 


AjA 

The  trend  which  is  exhibited  in  Case  2  also  appears  in  Cases  3,  4, 
and  5.  For  example,  at  r/r  ■  2,  the  peak  stress  in  ait  cases  lies  inside 
the  yield  surface.  Frequently,  at  greater  ranges  only  the  peak  or  only  the 
toe  of  the  stress  wave  lies  on  the  yield  surface.  A  Major  feature  of  the 
present  calculations  is  the  absence  of  sustained  plastic  deformation  on  loading. 

In  contrast,  unloading  is  frequently  accompanied  by  significant 

plastic  deformation.  The  outward  displacement  which  accompanies  any  compression 

wave  caused  o.  to  unload  more  rapidly  than  a  .  The  result  is  that  the 
o  r 

stress  point  in  a  /o.  plane  tends  away  from  the  hydros  tat  on  the  c  side 
until  It  encounters  the  yield  surface.  During  subsequent  plastic  unloading, 

0r  remains  the  major  principal  stress.  This  behavior  is  distinctly  different 
from  the  one~dimensiona1  strain  case,  in  which  the  major  principal  stress 
during  plastic  unloading  is  at  some  stages  perpendicular  to  the  direction  of 
wav_  propagation. 

PI LATENCY  AND  t/v 

Plastic  deformation  in  Case  4  is  associated  with  di latency,  or  plastic 
volumetric  expansion.  By  far  the  largest  amount  of  plastic  deformation,  and 
hence  di latency,  takes  place  during  unloading.  This  is  illustrated  in  Fig¬ 
ures  5*9  and  5*10  where  the  P/p  curves  for  Case  4  are  compared  with  those  In 

Case  2  at  r/r  ■  2  and  3.  Dilatency  causes  the  unloading  P/p  curve  to  lie  above 
o 

the  loading  curve.  Thus,  during  most  of  the  unloading  phase,  the  bulk  modulus 
is  less  than  either  the  loading  or  the  unloading  bulk  moduli  in  the  absence  of 
plasticity.  This  behavior,  in  conjunction  with  Inelastic  deformation  in  shear, 
reduces  the  speed  of  the  unloading  wave  and  therefore  contributes  to  a  slaver 
rate  of  attenuation  of  the  peak  or  for  Case  4  than  for  other  inelastic  cases. 

Although  the  unloading  P/p  curve  is  above  the  loading  curve,  there 

is  net  energy  dissipation  in  the  material.  This  is  demonstrated  in  Figures  5-9 

and  5*10,  where  the  areas  under  the  shear  stress/strain  curves  at  r/r  »  2 

o 

and  3,  representing  absorption  of  inelastic  energy,  are  greater  than  the  areas 
between  the  loading  and  unloading  P/p  curves. 
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a.  MEAN  STRESS/STRAIN  RELATION  AJAIOt) 


0  -0.01  -0.02  -0.03  -0.04 


b.  DEVIATOR  1C  STRESS/STRAIN  RELATION 


FIGURE  5-10.  ENERGY  DISSIPATION  AT  r/r  -  3 

o 
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GROWTH  OF  CAVITY 
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The  amount  of  cavity  expansion  is  a  reflection  of  the  material 
properties  and  the  amount  of  energy  originally  available  to  propagate  the 
wave.  The  displacement  of  the  cavity  surface  in  Cases  1-5  is  shown  as  a 
function  of  time  in  Figure  5*11.  Plastic  deformation  in  Cases  3*5  is  appar¬ 
ently  responsible  for  the  large  expansion.  One  reason  that  the  expansion 
in  Case  3  is  less  than  in  Cases  4  and  5  is  that  the  yield  criterion  in  Case  3 
permits  hoop  stresses  to  tend  toward  tension  whereas  the  criteria  in  Cases  A 
and  5  maintain  compressive  hoop  stress  at  the  cavity  wall.  The  second  reason 
is  that  mere  yielding  occurs  in  Cases  4  and  5  because  the  yield  criteria  are 
satisfied  at  a  lower  stress  level. 


FIGURE  5-11.  DISPLACEMENT  OF  CAVITY  WALL 
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ATTENUATION  RATES 

The  rate  at  which  peak  radial  stress  attenuates  with  range  is  a 
common  index  of  the  material  properties  and  wave  behavior.  The  stress  at 
the  front  of  a  spherically-diverging  step  wave  propagating  in  an  elastic 
material  attenuates  at  a  rate  proportional  to  r  ‘  (Reference  16).  inelas¬ 
ticity  may  accelerate  attenuation  by  allowing  unloading  signals  to  overtake 
the  front  and  by  absorbing  energy  of  the  wave. 

The  interpretation  of  attenuation  rates  in  the  present  work  is 
complicated  since  each  case  has  a  different  amount  of  energy  initially 
available  to  propagate  the  wave.  This  happens  because  a  common  pressure/time 
history  acts  through  different  displacements  of  the  cavity  wall.  The  amount 
of  energy  initially  available  to  the  wave  in  each  case  is  shown  in  Table  5-1. 


TABLE  5-1.  ENERGY  AVAILABLE  TO 
PROPAGATE  WAVE 


Case 

u/uo 

1 

1 

2 

1.03 

3 

1.52 

k 

2.61 

5 

2.69 

6 

8.95 

7 

**.25 

Uq  -  0.127  x  106  in. -lb 


The  attenuation  rates  of  peak  radial  stress  in  Cases  1-5  are  shown 
in  Figure  5-12.  Comparison  between  Cases  1  and  2  indicates  that  adding 
volumetric  hysteresis  strongly  increases  the  attenuation  rate.  This  is 
because  volumetric  hysteresis  absorbs  some  of  the  energy  of  the  wave  and 
because  the  unloading  signals  travel  faster  than  the  loading  signals,  over¬ 
taking  the  peak  and  degrading  it.  Interpretation  of  Case  is  complicated 
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because  there  is  considerably  more  energy  initially  available  to  propagate 
the  wave  than  in  Cases  1  and  2,  but  the  material  of  Case  4  is  tlso  much  more 
dissipative.  The  slow  attenuation  rate  of  the  peak  in  Case  1  appears  not 
to  be  due  to  the  extra  energy  in  the  wave,  which  is  primarily  associated  with 
the  portion  of  the  wave  behind  the  peak,  but  instead  to  the  reduction  in 
the  speed  of  unloading  signals  caused  by  dilstency.  This  is  due  to  reduction 
in  the  effective  bulk  modulus  caused  by  di latency  as  illustrated  in  Fig¬ 
ures  5*9a  and  5"10a. 


FIGURE  5-12.  EFFECT  OF  MATERIAL  PROPERTIES  AND 
IMPULSE  ON  ATTENUATION  RATE 
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Further  support  for  the  hypothesis  that  the  speed  of  unloading  waves 
dominates  attenuation  rates  comes  from  comparing  the  attenuation  rates  in 
Cases  2,  3,  and  5  with  that  in  Case  4.  This  represents  a  correlation  between 
unloading  wave  speeds  and  attenuation  rates.  There  is  no  such  simple  cor¬ 
relation  between  amount  of  energy  initially  in  the  wave  and  attenuation  rate, 
as  comparison  among  Cases  2,  4,  and  5  shows. 

The  most  dramatic  illustration  of  the  effect  of  volumetric  properties 
on  attenuation  rate  is  given  by  Case  4,  which  differs  from  Case  5  only  in 
having  plastic  dilatation.  At  r/rQ  «  4,  or  in  Case  4  is  twice  as  great  as 
that  in  Case  5. 


EFFECT  OF  PULSE  SHAPE  ON  WAVE  PROPERTIES 

To  test  the  hypothesis  that  the  rates  of  attenuation  depend  heavily 
on  the  unloading  waves  catching  up  with  the  peak  and  to  widen  the  scope  of  the 
study  in  general,  two  calculations  were  performed  using  input  pulses  of  longer 
duration.  The  input  pulse  used  in  Cases  6  and  7  is  illustrated  in  Figure  5-1b. 
The  peak  stress  on  the  cavity  wal i  is  the  same  as  in  Cases  1  through  5,  but  the 
peak  has  a  duration  of  7  x  10^  sec.  The  shape  of  the  unloading  pulse  is  similar 
to  Cases  1  through  5,  and  the  impulse  is  about  two  times  greater.  The  material 
properties  assumed  in  Case  6  are  the  same  as  those  in  Case  4  while  those  in 
Case  7  are  the  same  as  those  in  Case  2. 

The  stress/time  histories  in  Figure  5~5  and  the  rate  of  attenu¬ 
ation  plot,  Figure  5“ 1 3 »  confirm  that  lengthening  the  duration  of  the  peak 
delays  degradation  of  the  pulse.  Interpretation  of  Case  7  is  difficult  because 
the  peak  is  associated  with  a  very  sharp  rise  and  decay,  which  is  the  least 
favorable  wave  shape  for  the  present  method  of  calculation.  Comparison  of 
Case  2  with  Case  7  suggests  that  the  unloading  wave  catches  up  with  the  peak  at 

about  r/r  -  3-4  in  Case  7  and  r/r  -  1-2  in  Case  2.  This  is  confirmed 

o  o 

by  considering  the  instant  in  time  at  which  the  unloading  wave  catches  up  with 
the  front  of  the  loading  wave.  The  maximum  speed  of  unloading  waves  is 


1.94  x  10^  in. /sec 
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Unloading  Is  Initiated  at  r/r  , 

at  r/r  »  Z,  at  t  7  Ct  °-6  *  -t  1  “  7*67  x  I0"6s«c  ThJ 

„M  .  0  '  *  '  7’6?  *  10  6  sec  ♦  3  |„  ,„5  ™s  sl9nal  Arrives 

M  »PP™xiM,.ly  corrects  n-/’-M  «  '»5  I  a. /sec  -  2J.2  x 

«-**«.  h  Case  2  .he  un  I”  *'  °f  «  -  sea. 

■  mediately.  «v«  overtakes  the  peek  almost 


FIGURE  5-13  ccccrr 
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IMPLICATIONS  FOR  MATERIAL  PROPERTY  TESTING 


A  mathematical  model  transforms  experimental  data  from  the  states  of 
stress  under  which  the  data  were  measured  to  states  which  occur  in  calculations. 
Since  these  states  may  differ  widely,  great  reliance  must  be  placed  on  the 
mathematical  model  to  do  the  transformation  correctly.  Much  of  the  doubt  as 
to  the  physical  meaning  of  calculations  is  due  to  uncertainty  about  the  ability 
of  current  models  to  do  this. 


The  risk  involved  in  using  a  mathematical  model  in  the  present  cal¬ 
culations  could  be  lessened  by  basing  the  form  of  the  model  and  the  specific 
material  property  coefficients  on  laboratory  experiments  in  which  the  states 
of  stress  are  close  to  those  encountered  in  the  calculation.  One  use  of  the 
present  calculations  is  to  guide  the  planning  of  laboratory  experiments  on  - 
plane  or  cylindrical  samples  which  support  spherical  wave  field  tests.  The 
present  calculations  indicate  the  following  two  important  areas  of  concentration. 

a.  P/u  relation  for  hydrostatic  loading  and  unloading  to  indicate 
how  the  amount  of  hysteresis  varies  with  y 

max 

b.  Stress/strain  relations  for  loading  programs  in  which  the  yield 
criterion  and  inelastic  deformation  investigated  on  unloading. 

See  Figure  6-1. 


FIGURE  6-1.  PROPOSED  LOADING 
PROGRAMS  FOR  INVESTIGATING 
YIELD  CRITERIA  AND  INELASTIC 
DEFORMATION  ON  UNLOADING 
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The  purpose  of  investigating  the  P/p  relation  under  hydrostatic 
loading  is  to  establish  the  amount  of  hysteresis  in  the  absence  of  dilatancy, 
which  is  known  to  develop  near  the  maximum  deviatoric  stress.  Having  established 
this,  the  amount  of  dilatancy  which  occurs  during  inelastic  deformation  on 
unloading  may  be  investigated.  The  results  will  shed  light  on  which  of  the 
two  flow  rules  investigated  in  the  present  study  is  the  more  appropriate,  and 
whether  the  yield  criterion  is  the  same  on  unloading  as  on  loading. 

These  suggestions  are  intended  to  improve  the  effectiveness  of  the 
type  of  mathematical  model  used  in  the  present  calculation.  It  is  possible 
that  such  experiments  would  expose  the  present  model  as  being  inadequate  or 
misleading  for  application  to  spherical  wave  situation.  In  this  event,  it 
would  become  necessary  to  investigate  new  types  of  models. 
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CONCLUSIONS 

The  features  of  the  present  mathematics)  mode)  which  most  strongly 
affect  the  results  are 

a.  The  amount  of  permanent  volumetric  compaction  during  hydrostatic 
loading  and  unloading, 

b.  The  amount  of  dilatancy  accompanying  inelastic  deformation, 

c.  The  yield  criterion  during  unloading. 

Laboratory  testing  should  be  concentrated  in  these  areas  to  support  the 
development  of  a  model  for  use  in  spherical  wave  calculation. 
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An  investigation  was  made  of  the  effect  of  rise  time  on  the  pulse 
shape  and  attenuation  rate  of  spherically  diverging  waves.  The  purpose  of 
the  study  was  to  determine  whether  the  rise  time  of  the  prescribed  stress  or 
velocity  boundary  condition  signif icantiy  affects  the  calculated  stress  pulse 
at  ranges  between  1  and  5  times  the  cavity  radius.  For  the  rise  times  and 
material  properties  of  interest  to  the  present  finite  element  calculations, 
the  study  indicates  that  the  true  rise  time  should  be  represented  as  accurately 
as  possible.  The  initial  rise  time  of  about  0.67  x  10  ^  sec,  which  was 
measured  in  the  Physics  International  Company  experiments,  was  chosen  as 
being  representative  of  shock  waves  generated  by  detonating  chemical  explosive. 
The  finite  element  mesh  size  and  integration  time  step  were  selected  so  as 
to  represent  this  initial  rise  time  as  accurately  as  possible. 

In  the  first  part  of  the  study,  use  was  made  of  Jeffreys'  solution 
(Reference  i 7)  for  a  step  pressure  pulse  on  the  wail  of  a  cavity  whose  radius 
is  r  .  The  particle  velocity  at  a  range  r  is  given  by  the  following 
expression: 


1 

z 


■°-°2*>  (r  "  x)  s,n  6  +  2r  eos/F  al  exp(-B) 

Gr  L  J 


(A-1) 
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Radial  particle  velocity 

Cavity  radius 

Range 

Applied  radial  stress  on  the  cavity  wali 
biiitational  wave  speed 
Shear  modulus 
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The  solution  applies  subject  to  the  restriction  that 


Although  this  solution  applies  to  the  case  of  a  step  load  in  a 
cavity,  it  can  be  used  to  obtain  an  estimate  of  the  solution  for  a  ramp  load 
by  summing  up  the  effects  of  many  small  step  loads,  as  illustrated  in 
Figure  A-1.  Several  cases  were  studied  using  Jeffreys'  solution.  One  of 
these  was  compared  with  the  corresponding  finite  element  solution. 


Case 

Solution  Method 

t 

_R 

A1 

Jeffreys 

0. 

A2 

Jeffreys  Adapted 

0.0055  sec 

A3 

Fini te  Element 

0.0055  sec 

b.  APPftOJUMTICM  TO  RAMP  LOAD 

AJA1057 


FIGURE  A-1.  INPUT  FOR  JEFFREYS'  SOLUTION 
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The  following  geometrical,  material,  and  loadirg  paianeters  are  common  to  all 
three  cases. 


c  ■  1.465  x  10^  in. /sec 

P 

G  -  2.013  x  106  psl 

r  *  100  in. 

o 

°0  -  0.04775  psi 

The  velocity  time  history  for  each  case  is  obtained  at  r/r  -  2.44.  The 

o 

results  which  are  plotted  in  Figure  A-2  are  typical  of  the  dramatic  effect 
of  rise  time  on  the  peak  particle  velocity  and  pulse  shape. 


- finite  eiement  ys.oo»SEc 

- JEFFREYS'  SOLUTION  •  0.005$  SEC 

- JEFFREYS'  SOLUTION  ^  - 0. 


FIGURE  A-2.  VELOCITY/TIME  HISTORY  AT 
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A  second  study  was  initiated  to  study  the  effect  of  rise  time  on 
the  shape  and  attenuation  rate  of  a  wave  having  the  same  initial  shape  as 
the  one  used  in  the  present  finite  element  calculations.  This  shape  is 
represented  in  the  analysis  by  a  linear  rise  to  a  peak  pressure  followed 
by  an  exponential  decay.  Figure  A-3  illustrates  the  shape  of  the  input 
pulse  and  the  geometrical  coordinates  of  the  analysis. 


FIGURE  A-3.  COORDINATES  AND  PRESSURE/TIME 
HISTORY  FOR  RAMP  LOAD  WITH 
EXPONENTIAL  DECAY 


46 


R-681 3-777 


*5*0pr  _ 


AjA 

Consider  the  configuration  shown  in  Figure  A-3.  A  spherical  cavity 
of  radius  rQ  in  an  infinite  body  initially  at  rest  is  subjected  to  an 
interna]  pressure  p(t)  at  the  boundary  r  ■  r  .  Under  the  assumed  condition 
of  spherical  symmetry  and  infinitesimal  displacements,  the  response  at  a 
point  P  at  an  arbitrary  position  (r,  8,  is  to  be  determined.  Due  to 
spherical  symmetry,  the  components  of  displacement  in  the  8  and  direc¬ 
tions,  v  and  w,  respectively,  vanish.  The  remaining  response  parameters 
are  functions  of  r  and  t  only.  Since  only  infinitesimal  displacements 
are  considered,  no  distinction  need  be  made  between  deformed  and  undeformed 
posi tions. 


The  only  nontrivial  equation  of  motion  is  that  in  the  radial 
direction  (for  a  general  theory  of  elastic  wave  propagation;  see  Kolsky, 
Reference  18)  and  when  expressed  in  terms  of  the  dilatational  displacement 
potential  function  i|/  it  becomes  simply  the  one-dimensional  spherical  wave 
equation 


V2ij) 


=  UL+231 

"  3r2  r3r 


1  a2f 

2  2” 
c  3t 

P 


(A-2) 


where  ip  «  4«(r,t) 


!r  <  r  <  » 

0 

0  <  t  <  ® 


and  0^  -  \/(X  +  2p)/p  is  the  dilatational  wave  speed  and  X,  p  are  the 

Lamd  constants.  The  dilatation  potential  •  is  related  to  the  radial  displace¬ 
ment  u  by 


3r 


u  ■ 
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The  radial  stress  or  and  tangential  stresses  o^,  o  are  obtained  from 


<v  -  <x  ♦  2U)  |a ♦  2X  a.  .  (x  +  2p)L|+S.|1 


3r 


.  \  u  .  ,  3u  2 (X  +  \t)  3i|<  .  ,  3  i]/ 

°8  '  °*  •  2<A  *  ">  7*  X  57  ■  - ? — 57  *  *  72 

d  r 


All  shear  stresses  are  zero  as  a  result  of  symmetry. 


At  the  boundary  r  ■  r  , 


-P(t) 


(A-4) 


or 


(X  +  2u) 

Sr2  r  3r 

d  r 


-p(t) 


r  ■  r 


(A-5) 


and  with  r  «,  $  is  taken  to  be  zero  since  the  medium  far  enough  away  from 
the  cavity  is  undisturbed  for  finite  times. 

The  Laplace  transform  >f  <p(r,t)  is  denoted  by 


i(r,p) 


l ■- 


<l>(r,t)dt 


Applying  this  operation  to  the  differential  equation  (A-2)  provides 


14+  2  di.  k2^  .  0 

dr2  r  dr  d 


(A-6) 
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where  k^  ■  p/cp.  The  solution  to  Equation  A-6  Is,  accordingly, 

♦(r.p)  -  MeL.V  +  iM.V 


where  A(p),  B(p)  are  arbitrary  functions  of  p.  The  boundary  condition  at 
r  -*■  »  suggests  taking  B(p)  »  0  so  that 


*<r,p>  -  *Me'kPr 


(A- 7) 


Application  of  the  Laplace  transform  to  Equation  A-5  yields 


().  +  2p)  +  y-  4  -  -P(p) 

dr 


(A-8) 


where  P(p)  is  the  Laplece  transform  of  P(t).  Equations  A-7  and  A-8  together 
determine  the  arbitrary  function  A(p)  and  It  Is  found  that 


^(r.p) 


-  -p(r  -  r  ) / c 
.  lo  P  (p)e  °  p 

P»  (p  +  p1 ) (p  +  p2) 


(A-9) 


where  p^  *  2o(l  +  i 3) 
p2  «  2o(l  -  i e) 


2, 

c  /r  c 
sop 


■  Vw 


and  cg  ■  v/u/p  is  the  shear  wave  speed.  A  formal  solution  to  the  problem 
is  therefore  given  by  the  Inverse  Laplace  transform  of  Equat’on  A-9: 


'l'(r.t) 
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where  C  denotes  a  suitable  contour  in  the  complex  p-piane.  Equation  A-9 
can  also  be  written  as 


*(r,p) 


ro  P(p)  1 
pr  p2  -  P,  [p  +  P, 


1  1  -»(r  -  ro)/cp 

p  *  P2J 


(A— 11) 


which  is  a  more  convenient  form  for  the  inversion  process. 
Consider  the  exponential  puise  given  by 


^  ■  a r  [M(t)  - H(t  -  **.>] 

o  o 

(A— 1 2) 

1  -  t/t  N  -b  (t  -  At  )/t 

*  “rnr-  ZJ  v  [Hit  -  «g  -  hu  -  to>] 

where  H(t)  is  the  Heaviside  step  function 


H(t) 


|0  lft<« 

( 1  if  t  >  0 


and  N  is  any  integer.  The  pressure  puise  given  Ly  Equation  A-12  is  called 
an  N-term  exponential  pressure  puise.  It  can  be  shown  that  the  Lapiace  trans¬ 
form  of  Equation  A-12  is 


-At  p 
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_  1 

1 
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Substituting  Equation  A-l 3  into  Equation  A-11  and  carrying  out  the  inversion 
operation  leads  finally  to 


^r,t*  “  2c*8pr  'm 


7'm[^i['P't  +  v"jHt,)'['P'(T  A'0> 
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(A-14) 


where  t  =  t  -  (r  -  r  )/c 

o  p 

P1  =  2a  (1  +  i  8) 

2. 

a  "  cs/roCp 

6  "  \fT/ cs)2  -  1 

and  Im  denotes  the  imaginary  part  of  a  complex  function. 
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Expressions  for  u,  ar,  and  o0  can  be  obtained  from  Equations  A-3, 
A-4,  and  A-14.  The  radlai  velocity  u  is  obtained  by  differentiating  u 
with  respect  to  time  (the  lengthy  algebra  Involved  will  not  be  reproduced 

here).  Note  that  the  radial  displacement  u  has  terms  of  the  form  (r  /r) 

2  .  ° 
and  (r  /r)  ,  while  the  radial  velocity  u,  the  radial  stress  o  ,  and  the 
o  r_ 

tangential  stress  o  have  in  addition  terms  of  the  form  (r  /r)  .  The 

0  o 

response  has  the  general  character  of  a  highly  damped  oscillatory  motion 
about  the  static  value.  The  damping  time  constant  T  is  given  by 


T 


2a& 


where  a 

6 


cs/ro 

CP/CS 


Hence,  T  Is  of  the  order  of  the  half  transit  time  t.j.  =  a/ cg,  i.e.,  the  time 
for  a  shear  disturbance  to  travel  a  distance  equal  to  the  radius  of  the  cavity. 
Near  the*,  wave  front,  the  stresses  decay  like  (rQ/r)  with  increasing  radial 
distances.  Behind  the  front,  the  decay  rate  is  somewhat  higher  depending  on 
the  value  of  the  time  after  wave  front  arrival,  t,  compared  to  T.  The  limit, 
of  course,  is  (rQ/rr  for  the  static  case.  Therefore,  the  peak  response 
depends  not  only  on  the  peak  input  (as  in  all  linear  systems)  but  also  on  how 
soon  the  peak  is  attained.  In  other  words,  if  the  input  pulse  reaches  its 
peak  in  a  time  smaller  than  the  characteristic  damping  time  T,  the  response 
peaks  are  expected  to  scale  according  to  (r/rQ)  ^ •  If  the  input  peak  is 
not  reached  until  after  several  multiples  of  T,  however,  the  response  peaks 
are  expected  to  scale  according  to  (r/rQ)  •  This  observation  is  verified 
by  the  results  obtained  for  the  specific  examples  considered  below. 


Numerical  results  are  obtained  for  the  exponential  input  pulses 
defined  by  the  fnllowinn  parameters. 
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al 

m 

0.65 

a2 

m 

0.35 

bl 

m 

10 

b2 

m 

0.1 

Po 

m 

-475,000  psi 

*o 

« 

40  x  10'6  sec 

Case 


Bl 

0. 

B2 

0.6667  x  10"6 

P3 

2.0  x  10"6  sec 

B4 

5.0  x  10  6  sec 
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These  parameters  describe  pressure  time  histories  on  the  surface  of  the  cavity 
which  cover  the  range  observed  in  the  present  finite  element  calculations. 
Although  the  input  pulse  in  all  the  finite  element  calculations  has  an  initial 
rise  time  of  0.6667  x  10  ^  sec  (Case  B2) ,  the  rise  time  lengthens  to 
2.0  x  10  b  to  3.0  x  10  6  at  r/rQ  “  4  to  5.  This  lengthening  is  due  to 
numerical  distortion  of  the  pulse  which  evidently  cannot  be  entirely  elimi¬ 
nated  in  spite  of  using  small  mesh  size  and  integration  time  step.  Thus, 

Cases  B'  through  B4  bound  the  possible  effects  on  pulse  shape  and  attenuation 
rate  which  can  be  attributed  to  rise  time  in  the  present  elastic  finite 
element  calculations. 

The  medium  surrounding  the  cavity  has  the  same  properties  as  the 
medium  in  Case  1  of  the  finite  element  calculations.  These  are 

B  *  Bq  =  1 .205  x  10b  ps i 

G  -  1 .000  x  106  psi 

P  -  0.000238  lb-sec2/in.i* 


53 


R-6613-777 


AjA 

The  p-  and  s-wave  speeds  in  the  material  are 

c  ■  8590  ft/sec 
P 

cg  *  5^00  ft/sec 

The  cavity  is  assumed  to  have  an  initial  radius  r  ■  1  in. 

The  responses  in  the  four  cases  are  shown  at  r/rQ  =  1,  2,  3,  and  4 
in  Figures  A-4  and  A-5.  The  computer  program  which  was  used  to  calculate  the 
response  automatically  defines  the  arrival  time  at  each  station  to  be  zero. 

Hence,  the  response  at  all  tour  ranges  of  a  given  case  appear  to  start  at 
zero,  whereas  in  a  physical  problem  each  would  have  different  times  of  arrival. 

Except  in  Case  B4,  when  the  portion  of  the  of/t  history  prior  to 

the  peak  is  slightly  nonlinear,  there  is  no  noticeable  distortion  of  the  pulse 

due  to  the  various  rise  time  studied.  The  rates  of  attenuation,  which  are 

shown  in  Figure  A-6,  differ  appreciably  due  to  the  different  rise  times.  For 

example,  at  r/r  ■  4,  0  /u  in  Case  B1  is  1.45  times  that  in  Case  B4.  This 
o  r  o 

effect  is  considered  to  be  too  big  to  be  ignored,  and  hence,  care  was  used 
In  representing  rise  times  in  the  finite  element  calculations. 

An  estimate  cf  the  numerical  error  present  in  the  finite  element 
calculations  is  indicated  by  comparing  the  attenuation  rate  found  in  Case  1 
with  the  rates  of  attenuation  found  in  the  closed  form  solutions. 

The  formal  integral  solution  is  obtained  for  spherical  waves  propagating 

from  a  spherical  cavity  subjected  to  a  pressure  pulse  with  finite  rise  time 

and  exponential  decay.  The  rate  of  attenuation  of  the  peak  stress  when  it 

propagates  into  the  medium  depends  on  the  magnitude  of  the  rise  time  relative 

to  the  transit  time  r  /c  ,  where  r  is  the  radius  of  the  cavity  and  c 

o  s  o  s 

is  the  shear  wave  speed.  Fcr  zero  rise  time,  the  attenuation  rate  is  propor¬ 
tional  to  r  \  whereas  for  rise  times  much  greater  than  r  /c  the  attenuation 

■3  os 

rate  ir  proportional  to  r  . 


54 


FIGURE  A-A.  RADIAL-STRESS/TIME  HISTORIES  FOR  VARIOUS  RISE  TIMES 
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FIGURE  A-5.  RADIAL-VELOCITY/TIME  HISTORIES  FOR  VARIOUS  RISE  TIMES 
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FIGURE  A-6.  EFFECT  OF  RISE  TIME  ON  ATTENUATION  RATE 
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